/*
 * Copyright (c) 2003-2008 Fred Hutchinson Cancer Research Center
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package edu.bu.manateem.msg.heatmap;

import java.io.IOException;
import java.io.Writer;
import java.util.Arrays;
import java.util.Comparator;
import java.util.List;

import org.apache.commons.lang.math.FloatRange;

/**
 * User: mbellew
 * Date: Jun 2, 2004
 * Time: 1:57:21 PM
 */
public class Spectrum
{
	// private static Logger _log = Logger.getLogger(Spectrum.class);
	
	/*
	 * this is just a marker class to distinguish between raw and centroided
	 * peaks
	 */
	public static class CentroidedPeak
		extends Peak
	{
		public CentroidedPeak(float mz,float intensity)
		{
			super(mz,intensity);
		}
		
		public CentroidedPeak(int scan,float mz,float intensity)
		{
			super(scan,mz,intensity);
		}
	}
	public static class Peak
	{
		public int scan = -1;
		public float mz;
		public float intensity;
		public float background = 0F;
		public float median = 0F;
		public boolean excluded = false;
		
		public Peak()
		{
			// empty block
		}
		
		public Peak(float mz,float intensity)
		{
			this.mz = mz;
			this.intensity = intensity;
		}
		
		public Peak(int scan,float mz,float intensity)
		{
			this.scan = scan;
			this.mz = mz;
			this.intensity = intensity;
		}
		
		public Peak(Peak p)
		{
			scan = p.scan;
			mz = p.mz;
			intensity = p.intensity;
			background = p.background;
			median = p.median;
		}
		
		@Override
		public boolean equals(Object o)
		{
			if ((null == o) || (o.getClass() != this.getClass()))
			{
				return false;
			}
			
			Peak p = (Peak) o;
			return (p.scan == scan) && (p.mz == mz)
				&& (p.intensity == intensity);
		}
		
		public float getBackground()
		{
			return background;
		}
		
		public float getIntensity()
		{
			return intensity;
		}
		
		public float getMedian()
		{
			return median;
		}
		
		public float getMz()
		{
			return mz;
		}
		
		public int getScan()
		{
			return scan;
		}
		
		@Override
		public int hashCode()
		{
			return scan ^ Float.floatToIntBits(mz)
				^ Float.floatToIntBits(intensity);
		}
		
		public void setBackground(float background)
		{
			this.background = background;
		}
		
		public void setIntensity(float intensity)
		{
			this.intensity = intensity;
		}
		
		public void setMedian(float median)
		{
			this.median = median;
		}
		
		public void setMz(float mz)
		{
			this.mz = mz;
		}
		
		public void setScan(int scan)
		{
			this.scan = scan;
		}
		
		@Override
		public String toString()
		{
			return scan + "\t" + mz + "\t" + intensity;
		}
	}
	
	/**
	 *for micromass, doesn't work for bruker
	 public static float[][] ExpandZeroes(float[][] spectrum)
	 {
	 final float ZERO = 0F; 			// set to -2 for debugging
	 final float ZEROFILL = 0F;     	// set to -1 for debugging
	 final int notZero = 0;
	 final int oneZero = 1;
	 final int twoZero = 2;

	 float[] mzIN = spectrum[0];
	 float[] intIN = spectrum[1];
	 FloatArray mzOUT = new FloatArray();
	 FloatArray intOUT = new FloatArray();

	 int state = notZero;
	 float interval = Math.min(0.5F, mzIN[1] - mzIN[0]);
	 for (int i=0 ; i<mzIN.length ; i++)
	 {
	 float mz = mzIN[i];
	 float f = intIN[i];
	 if (f != 0)
	 {
	 state = notZero;
	 mzOUT.add(mz);
	 intOUT.add(f);
	 }
	 else if (state == notZero)
	 {
	 state = oneZero;
	 mzOUT.add(mz);
	 intOUT.add(ZERO);
	 }
	 else
	 {
	 state = twoZero;
	 if (i >= 2)
	 interval = mzIN[i-1] - mzIN[i-2];
	 float mzFrom = mzIN[i-1];
	 float mzTo = mzIN[i] - interval/2;
	 for (float m = mzFrom+interval ; m < mzTo ; m+=interval)
	 {
	 mzOUT.add(m);
	 intOUT.add(ZEROFILL);
	 }
	 mzOUT.add(mz);
	 intOUT.add(ZERO);
	 }
	 }
	 return new float[][] {mzOUT.toArray(null), intOUT.toArray(null)};
	 }
	 */
	
	public static final float HYDROGEN_ION_MASS =
		(float) (1.0078250 - 5.485e-4); // h - electron
	
	public static final double LN2 = Math.log(2.0);
	
	public static Comparator<Peak> comparePeakIntensityDesc =
		new Comparator<Peak>()
		{
			public int compare(Peak a,Peak b)
			{
				// if (a == null || b == null) return a==null ? 1 : -1;
				return Spectrum._compareDesc(a.intensity,b.intensity);
			}
		};
	
	public static Comparator<Peak> comparePeakMzAsc = new Comparator<Peak>()
	{
		public int compare(Peak a,Peak b)
		{
			return Spectrum._compareAsc(a.mz,b.mz);
		}
	};
	
	public static Comparator<Peak> comparePeakScanAsc = new Comparator<Peak>()
	{
		public int compare(Peak a,Peak b)
		{
			return Spectrum._compareAsc(a.scan,b.scan);
		}
	};
	
	static float[][] _poisson = new float[6400 / 10][6];
	
	static
	{
		for (int i = 0; i < Spectrum._poisson.length; i++)
		{
			int mass = i * 10 + 5;
			float mu = mass * 0.0005556F; // martin's empircal lambda = 1/1800
											// (was mass * 0.000489F)
			double muExp = Math.exp(-mu);
			Spectrum._poisson[i][0] = (float) (1 * muExp / 1);
			Spectrum._poisson[i][1] = (float) (mu * muExp / 1);
			Spectrum._poisson[i][2] = (float) (Math.pow(mu,2) * muExp / 2);
			Spectrum._poisson[i][3] = (float) (Math.pow(mu,3) * muExp / 6);
			Spectrum._poisson[i][4] = (float) (Math.pow(mu,4) * muExp / 24);
			Spectrum._poisson[i][5] = (float) (Math.pow(mu,5) * muExp / 120);
			float s =
				Spectrum._poisson[i][0] + Spectrum._poisson[i][1]
					+ Spectrum._poisson[i][2] + Spectrum._poisson[i][3]
					+ Spectrum._poisson[i][4] + Spectrum._poisson[i][5];
			Spectrum._poisson[i][0] /= s;
			Spectrum._poisson[i][1] /= s;
			Spectrum._poisson[i][2] /= s;
			Spectrum._poisson[i][3] /= s;
			Spectrum._poisson[i][4] /= s;
			Spectrum._poisson[i][5] /= s;
		}
	}
	
	// from article in Electrophoresis 1999 by S. Gay et al
	static double[] m0 = new double[]
		{
			0.0576,-0.2553,0.5827,-0.8436,0.6182
		};
	
	// static class FloatImage extends PlanarImage
	// {
	// WritableRaster raster;
	//
	//
	// FloatImage(int w, int h, DataBuffer buf)
	// {
	// tileWidth = width = w;
	// tileHeight = height = h;
	// tileGridXOffset = tileGridYOffset = minX = minY = 0;
	// assert buf.getSize() == w * h;
	// sampleModel =
	// RasterFactory.createBandedSampleModel(DataBuffer.TYPE_FLOAT, w, h, 1);
	// raster = RasterFactory.createWritableRaster(sampleModel, buf, new
	// Point(0, 0));
	// }
	//
	//
	// public Raster getData()
	// {
	// return raster;
	// }
	//
	//
	// public Raster getTile(int tileX, int tileY)
	// {
	// return raster;
	// }
	// }
	
	// public static float[] FFTsmooth(float[] x, double smoothfactor, boolean
	// cliff)
	// {
	// try
	// {
	// DataBuffer buf = new DataBufferFloat(x, x.length);
	// FloatImage img = new FloatImage(x.length, 1, buf);
	// RenderedOp opDFT = JAI.create("DFT", new ParameterBlock()
	// .add(DFTDescriptor.SCALING_UNITARY)
	// .add(DFTDescriptor.REAL_TO_COMPLEX)
	// .addSource(img));
	//
	// float[] fft = ((DataBufferFloat)
	// opDFT.getData().getDataBuffer()).getData();
	// int complexLen = fft.length / 2;
	//
	// // UNDONE cache, this seems overexpensive
	// /*
	// int len = complexLen / 2;
	// float[] g = new float[len];
	// double s = len / smoothfactor;
	// double s2 = -1.0 / (2.0 * s * s);
	// int z = cliff ? (int) s : 0;
	// for (int i = 0; i < z; i++)
	// g[i] = 1F;
	// for (int i = 0; i < len - z; i++)
	// {
	// double d = Math.exp((double)i * i * s2);
	// if (d < Float.MIN_VALUE) // arbitrary cut-off for small..
	// break;
	// g[z + i] = (float) d;
	// }
	// */
	// int len = complexLen / 2;
	// float[] g = NormalP(len, 1/smoothfactor, false, false, null);
	//
	// int off = 2;
	// for (int i = 0; i < len; i++)
	// {
	// fft[off++] *= g[i]; // Real
	// fft[off++] *= g[i]; // Complex
	// }
	// for (int i = len - 2; i >= 0; i--)
	// {
	// fft[off++] *= g[i]; // Real
	// fft[off++] *= g[i]; // Complex
	// }
	// assert off == fft.length;
	//
	// RenderedOp opIDFT = JAI.create("IDFT", new ParameterBlock()
	// .add(DFTDescriptor.SCALING_UNITARY)
	// .add(DFTDescriptor.COMPLEX_TO_REAL)
	// .addSource(opDFT));
	// DataBufferFloat bufReturn = (DataBufferFloat)
	// opIDFT.getData().getDataBuffer();
	// return bufReturn.getData();
	// }
	// catch (Exception ex)
	// {
	// ex.printStackTrace();
	// return null;
	// }
	// }
	
	static double[] m1 = new double[]
		{
			-0.1207,0.4688,-0.7373,0.3845,0.2701
		};
	
	/*
	 * Computes points in Gaussian distribution function, over range [0,1]
	 * with mean=0.5
	 * @param len number of points to compute in the distrubition
	 * @param sigma sqrt(variance)
	 * @param g optional array to copy result into
	 * @return
	 * private static float[] NormalD(int len, double sigma, float[] g)
	 * {
	 * assert null == "NYI";
	 * return null;
	 * }
	 */

	static double[] m2 = new double[]
		{
			0.0569,-0.1173,-0.0838,0.3068,0.0865
		};
	
	static double[] m3 = new double[]
		{
			0.0296,-0.1293,0.1330,0.1151,0.0204
		};
	
	static double[] m4 = new double[]
		{
			-0.0101,-0.0118,0.0789,0.0292,0.0040
		};
	
	static double[] m5 = new double[]
		{
			-0.0132,0.0448,0.0256,0.0079,0.0008
		};
	
	static float[][] _gay = new float[6400 / 10][6];
	
	static
	{
		for (int i = 0; i < Spectrum._gay.length; i++)
		{
			int m = i * 10 + 5;
			double x = (m - 800.0) / 2200.0;
			double x2 = x * x;
			double x3 = x2 * x;
			double x4 = x3 * x;
			Spectrum._gay[i][0] =
				(float) (Spectrum.m0[0] * x4 + Spectrum.m0[1] * x3
					+ Spectrum.m0[2] * x2 + Spectrum.m0[3] * x + Spectrum.m0[4]);
			Spectrum._gay[i][1] =
				(float) (Spectrum.m1[0] * x4 + Spectrum.m1[1] * x3
					+ Spectrum.m1[2] * x2 + Spectrum.m1[3] * x + Spectrum.m1[4]);
			Spectrum._gay[i][2] =
				(float) (Spectrum.m2[0] * x4 + Spectrum.m2[1] * x3
					+ Spectrum.m2[2] * x2 + Spectrum.m2[3] * x + Spectrum.m2[4]);
			Spectrum._gay[i][3] =
				(float) (Spectrum.m3[0] * x4 + Spectrum.m3[1] * x3
					+ Spectrum.m3[2] * x2 + Spectrum.m3[3] * x + Spectrum.m3[4]);
			Spectrum._gay[i][4] =
				(float) (Spectrum.m4[0] * x4 + Spectrum.m4[1] * x3
					+ Spectrum.m4[2] * x2 + Spectrum.m4[3] * x + Spectrum.m4[4]);
			Spectrum._gay[i][5] =
				(float) (Spectrum.m5[0] * x4 + Spectrum.m5[1] * x3
					+ Spectrum.m5[2] * x2 + Spectrum.m5[3] * x + Spectrum.m5[4]);
			float s =
				Spectrum._gay[i][0] + Spectrum._gay[i][1] + Spectrum._gay[i][2]
					+ Spectrum._gay[i][3] + Spectrum._gay[i][4]
					+ Spectrum._gay[i][5];
			Spectrum._gay[i][0] /= s;
			Spectrum._gay[i][1] /= s;
			Spectrum._gay[i][2] /= s;
			Spectrum._gay[i][3] /= s;
			Spectrum._gay[i][4] /= s;
			Spectrum._gay[i][5] /= s;
		}
	}
	
	static int _compareAsc(float a,float b)
	{
		return a == b
			? 0
			: a < b
				? -1
				: 1;
	}
	
	static int _compareDesc(float a,float b)
	{
		return a == b
			? 0
			: a < b
				? 1
				: -1;
	}
	
	static float[] _gay(float m)
	{
		int i = ((int) m - 5) / 10;
		i = Math.max(0,Math.min(Spectrum._gay.length - 1,i));
		return Spectrum._gay[i];
	}
	
	/**
	 * Return the indexes of the peaks in descending order of theoretical intensity.
	 * Relies on no exact ties from Poisson();
	 * @param mass
	 * @return
	 */
	public static int[] calcIdealPeakIntensityOrderDesc(float mass)
	{
		// Compare intensities of the peak that has theoretical max intensity
		float[] p = Spectrum.Poisson(mass);
		float[] pSorted = new float[p.length];
		System.arraycopy(p,0,pSorted,0,p.length);
		Arrays.sort(pSorted);
		// sorted in ASCENDING ORDER
		
		int[] result = new int[p.length];
		for (int i = 0; i < p.length; i++)
		{
			float match = pSorted[pSorted.length - i - 1];
			for (int j = 0; j < p.length; j++)
			{
				if (p[j] == match)
				{
					result[i] = j;
					break;
				}
			}
		}
		return result;
	}
	
	/**
	 * Return the index of the peak that, according to the Poisson distribution for this mass, should be most intense
	 * @param mass
	 * @return
	 */
	public static int calcMaxIdealPeakIndex(float mass)
	{
		// Compare intensities of the peak that has theoretical max intensity
		int theoreticalMaxIndex = 0;
		float[] p = Spectrum.Poisson(mass);
		float maxTheoreticalInt = p[0];
		for (int i = 1; i < p.length; i++)
		{
			if (p[i] > maxTheoreticalInt)
			{
				maxTheoreticalInt = p[i];
				theoreticalMaxIndex = i;
			}
		}
		return theoreticalMaxIndex;
	}
	
	public static double ChooseThreshold(double[] x,double f)
	{
		{
			double m = Spectrum.MedianSampled(x,true);
			return m * f;
		}
	}
	
	public static float[][] CombineRawSpectra(float[] template,
		List<float[][]> spectra,FloatRange rangeIN)
	{
		int from = Spectrum.find(template,rangeIN.getMinimumFloat());
		if ((from > 0) && (template[from] > rangeIN.getMinimumFloat()))
		{
			from--;
		}
		int to = Spectrum.find(template,rangeIN.getMaximumFloat());
		if (to == template.length)
		{
			to--;
		}
		float[] mzArray = new float[to - from + 1];
		System.arraycopy(template,from,mzArray,0,mzArray.length);
		
		// now combine intensity values
		float[] iArray = new float[mzArray.length];
		for (int s = 0; s < spectra.size(); s++)
		{
			float[][] spectrum = spectra.get(s);
			from = Spectrum.find(spectrum[0],rangeIN.getMinimumFloat());
			to = Spectrum.find(spectrum[0],rangeIN.getMaximumFloat());
			if (spectrum[0][to] > rangeIN.getMaximumFloat())
			{
				--to;
			}
			for (int i = from,j = 0; i < to; i++)
			{
				float f = spectrum[0][i];
				while (mzArray[j] < f)
				{
					j++;
				}
				// which is closer mzArray[j] or mzArray[j-1];
				if ((j == 0) || (mzArray[j] - f < f - mzArray[j - 1]))
				{
					iArray[j] += spectrum[1][i];
				}
				else
				{
					iArray[j - 1] += spectrum[1][i];
				}
			}
		}
		int sz = spectra.size();
		for (int i = 0; i < iArray.length; i++)
		{
			iArray[i] /= sz;
		}
		return new float[][]
			{
				mzArray,iArray
			};
	}
	
	public static void CopyToTSV(float[][] spectrum,Writer out,boolean useHeader)
		throws IOException
	{
		if (useHeader)
		{
			out.write("mz\tintensity\n");
		}
		for (int i = 0; i < spectrum[0].length; i++)
		{
			out.write(String.valueOf(spectrum[0][i]));
			out.write('\t');
			out.write(String.valueOf(spectrum[1][i]));
			out.write('\n');
		}
	}
	
	public static double Correlation(float[] X,float[] Y)
	{
		assert X.length == Y.length;
		if (X.length != Y.length)
		{
			throw new IllegalArgumentException(
				"Spectrum.Correlation(), arrays need be equal length");
		}
		int n = X.length;
		double Sx = 0.0,Sx2 = 0.0,Sy = 0.0,Sy2 = 0.0,Sxy = 0.0;
		for (int i = 0; i < n; i++)
		{
			double x = X[i],y = Y[i];
			Sx += x;
			Sx2 += x * x;
			Sy += y;
			Sy2 += y * y;
			Sxy += x * y;
		}
		return (n * Sxy - Sx * Sy)
			/ Math.sqrt((n * Sx2 - Sx * Sx) * (n * Sy2 - Sy * Sy));
	}
	
	static int find(float[] s,float f)
	{
		int i = Arrays.binarySearch(s,f);
		i = i < 0
			? -(i + 1)
			: i;
		i = Math.min(s.length - 1,i);
		return i;
	}
	
	/**
	 * For both micromass and bruker, intervals between readings are proportional to sqrt(mz).
	 * This function infers the sampling frequency from a spectrum.
	 *
	 * r.min should be an actual start value from a real spectrum (e.g. not just 400)
	 *
	 * This is expensive so it's done once per MSRun, see MSRun.getTemplateSpectrum().
	 * Since we do it once, don't try to optimize too much.
	 */
	public static float[] GenerateSpectrumTemplate(float[] mzArray,FloatRange r)
	{
		// convert to sqrt(mz) domain
		double start_sqrt = Math.sqrt(r.getMinimumDouble());
		double end_sqrt = Math.sqrt(r.getMaximumDouble());
		double[] mzArray_sqrt = new double[mzArray.length];
		for (int i = 0; i < mzArray.length; i++)
		{
			mzArray_sqrt[i] = Math.sqrt(mzArray[i]);
		}
		
		/*
		 * double gaps[] = new double[mzArray.length-1];
		 * int lenGaps = 0;
		 * for (int i=1 ; i<mzArray.length ; i++)
		 * {
		 * double d = mzArray_sqrt[i] - mzArray_sqrt[i-1];
		 * if (d > 0 && d < .025)
		 * gaps[lenGaps++] = d;
		 * }
		 * Arrays.sort(gaps, 0, lenGaps);
		 */
		// find smallest gap
		double min = 0.001;
		for (int i = 1; i < mzArray_sqrt.length; i++)
		{
			double d = mzArray_sqrt[i] - mzArray_sqrt[i - 1];
			if (d < min)
			{
				min = d;
			}
		}
		
		// find average of small gaps
		double sum = 0.0;
		int count = 0;
		double cutoff = min * 1.5;
		for (int i = 1; i < mzArray_sqrt.length; i++)
		{
			double d = mzArray_sqrt[i] - mzArray_sqrt[i - 1];
			if (d < cutoff)
			{
				sum += d;
				count++;
			}
		}
		
		double interval_sqrt = sum / count;
		// double interval_upper = interval_sqrt * 1.5;
		double interval_lower = interval_sqrt * 0.5;
		
		// We have a really good interval value, but we'll drift from reported
		// mz values if we use it directly.
		// "merge" input data and calculate missing values
		int lenTemp =
			(int) Math.floor((end_sqrt - start_sqrt) / interval_sqrt) + 200;
		double[] template_sqrt = new double[lenTemp];
		
		int dst = 0,src = 0;
		mzArray_sqrt[0] = start_sqrt;
		mzArray_sqrt[mzArray_sqrt.length - 1] = end_sqrt;
		double prev = start_sqrt;
		for (; src < mzArray_sqrt.length; src++)
		{
			double next = mzArray_sqrt[src];
			double gap = next - prev;
			int c = (int) ((gap + interval_lower) / interval_sqrt);
			assert (c < 2)
				|| ((gap / c > interval_sqrt * 0.9) && (gap / c < interval_sqrt * 1.1));
			for (int i = 1; i < c; i++)
			{
				template_sqrt[dst++] = prev + i * gap / c;
			}
			template_sqrt[dst++] = prev = next;
		}
		int len = dst;
		
		// convert back to mz domain
		float[] template = new float[len];
		for (int i = 0; i < len; i++)
		{
			double mz = template_sqrt[i] * template_sqrt[i];
			template[i] = (float) (Math.round(10000.0 * mz) / 10000.0);
			// System.out.println(template[i]);
		}
		
		return template;
	}
	
	public static float[] getRow(float[][] m,int r,float[] out)
	{
		out = Spectrum.realloc(out,m.length);
		for (int s = 0; s < m.length; s++)
		{
			out[s] = m[s][r];
		}
		return out;
	}
	
	public static double HellingerDistance(float[] p,float[] q)
	{
		assert p.length == q.length;
		double d = 0.0;
		for (int i = 0; i < p.length; i++)
		{
			double x = Math.sqrt(p[i]) - Math.sqrt(q[i]);
			d += x * x;
		}
		return d;
	}
	
	private static float[] Interpolate(float[] m,int len,int windowSize,
		float[] result)
	{
		windowSize = Math.min(len,windowSize);
		try
		{
			result = Spectrum.realloc(result,len);
			// fill non-interpolated part
			Arrays.fill(result,0,windowSize,m[0]);
			Arrays.fill(result,result.length - windowSize,result.length,
				m[m.length - 1]);
			// fill interpolated part
			for (int w = 0; w < m.length - 1; w++)
			{
				int windowStart = len * w / m.length + windowSize / 2;
				int windowEnd = len * (w + 1) / m.length + windowSize / 2;
				int windowLen = windowEnd - windowStart;
				double v = m[w];
				double d = (m[w + 1] - m[w]) / windowLen;
				for (int i = 0; i < windowLen; i++)
				{
					result[windowStart + i] = (float) (v + d * i);
				}
			}
			return result;
		}
		catch (RuntimeException x)
		{
			x.printStackTrace();
			throw x; // just a place to put a breakpoint
		}
	}
	
	public static double KLDistanceSymmetric(float[] signal1,int off1,
		float[] signal2,int off2,int len)
	{
		double sum1 = 0.0,sum2 = 0.0;
		for (int i = 0; i < len; i++)
		{
			sum1 += signal1[off1 + i];
			sum2 += signal2[off2 + i];
		}
		
		double diff = 0.0;
		double min = 0.001 / len;
		for (int i = 0; i < len; i++)
		{
			double s1 = Math.max(signal1[off1 + i] / sum1,min);
			double s2 = Math.max(signal2[off2 + i] / sum2,min);
			diff += s1 * Math.log(s1 / s2) + s2 * Math.log(s2 / s1);
		}
		return diff / Math.log(2.0);
	}
	
	public static float KLGayDistance(float m,float[] signal)
	{
		float[] p = Spectrum._gay(m);
		float diff = 0.0F;
		for (int k = 0; k < 4; k++)
		{
			float f0 = p[k];
			float f1 = signal[k];
			diff += f0 * Math.log(f0 / f1) + f1 * Math.log(f1 / f0);
		}
		return diff;
	}
	
	public static float KLPoissonDistance(float mass,float[] q)
	{
		float[] p = Spectrum.Poisson(mass);
		assert q.length == p.length;
		
		double diff = 0.0;
		double sumP = 0.0;
		double sumQ = 0.0;
		int pqMinLength = Math.min(p.length,q.length);
		for (int k = 0; k < pqMinLength; k++)
		{
			diff += p[k] * Math.log((double) p[k] / q[k]);
			sumP += p[k];
			sumQ += q[k];
		}
		double kl = diff / Spectrum.LN2;
		
		assert kl > -0.0001;
		assert (Math.abs(sumP - 1.0) < 0.001) && (Math.abs(sumQ - 1.0) < 0.001);
		
		kl = Math.max(kl,0);
		return (float) (diff / Math.log(2.0));
	}
	
	public static float KLPoissonDistanceSymmetric(float mass,float[] signal)
	{
		float[] p = Spectrum.Poisson(mass);
		float diff = 0.0F;
		for (int k = 0; k < 4; k++)
		{
			float f0 = p[k];
			float f1 = signal[k];
			diff += f0 * Math.log(f0 / f1) + f1 * Math.log(f1 / f0);
		}
		return (float) (diff / Math.log(2.0));
	}
	
	public static void main(String[] args)
	{
		float[] g10 = Spectrum.NormalP(11,.01,false,false,null);
		float[] g20 = Spectrum.NormalP(21,.01,false,false,null);
		for (int i = 0; i < g10.length; i++)
		{
			System.out.println("" + (((double) i) / (g10.length - 1)) + "\t"
				+ g10[i]);
		}
		System.out.println();
		for (int i = 0; i < g20.length; i++)
		{
			System.out.println("" + (((double) i) / (g20.length - 1)) + "\t"
				+ g20[i]);
		}
		
		assert 2F == Spectrum.Median(1,2,3);
		assert 2F == Spectrum.Median(1,3,2);
		assert 2F == Spectrum.Median(2,1,2);
		assert 2F == Spectrum.Median(2,1,3);
		assert 2F == Spectrum.Median(2,3,1);
		assert 2F == Spectrum.Median(3,2,1);
	}
	
	/*
	 * public static boolean IsLockSpray(float[][] fullSpectrum)
	 * {
	 * float[][] spectrum = ResampleSpectrum(fullSpectrum, new
	 * FloatRange(780.0F, 7900.F), 36, false);
	 * spectrum[1] = Spectrum.FFTsmooth(spectrum[1], 8, false);
	 * Peak[] peaks = Spectrum.PickPeaks(spectrum, 20.0F);
	 * Arrays.sort(peaks, Spectrum.comparePeakIntensityDesc);
	 * if (peaks.length < 3)
	 * return false;
	 * if (Math.abs(peaks[0].mz - 785.8333) > .06)
	 * return false;
	 * if (Math.abs(peaks[1].mz - 786.3333) > .06)
	 * return false;
	 * if (Math.abs(peaks[2].mz - 786.8333) > .06)
	 * return false;
	 * return true;
	 * }
	 */

	public static double Median(double[] x,int start,int len,boolean fABS,
		double[] t)
	{
		assert (null == t) || (t.length >= len);
		t = Spectrum.realloc(t,len);
		if (!fABS)
		{
			System.arraycopy(x,start,t,0,len);
		}
		else
		{
			for (int i = 0; i < len; i++)
			{
				t[i] = Math.abs(x[i + start]);
			}
		}
		Arrays.sort(t);
		return t[len / 2];
	}
	
	public static float Median(float a,float b,float c)
	{
		if (a < b)
		{
			return c < a
				? a
				: b < c
					? b
					: c;
		}
		else
		{
			return a < c
				? a
				: c < b
					? b
					: c;
		}
	}
	
	public static float Median(float[] x,int start,int len,boolean fABS,
		float[] t)
	{
		assert (null == t) || (t.length >= len);
		t = Spectrum.realloc(t,len);
		if (!fABS)
		{
			System.arraycopy(x,start,t,0,len);
		}
		else
		{
			for (int i = 0; i < len; i++)
			{
				t[i] = Math.abs(x[i + start]);
			}
		}
		Arrays.sort(t);
		return t[len / 2];
	}
	
	public static double MedianSampled(double[] x,boolean fABS)
	{
		double p = 0.5;
		int len = x.length;
		int sample = 100 + (int) Math.ceil(Math.sqrt(len));
		
		double[] copy = null;
		if (sample < len / 4)
		{
			copy = new double[sample];
			for (int i = 0; i < sample; i++)
			{
				int r = (int) Math.floor(Math.random() * len);
				copy[i] = fABS
					? Math.abs(x[r])
					: x[r];
			}
		}
		else
		{
			copy = new double[len];
			for (int i = 0; i < len; i++)
			{
				copy[i] = fABS
					? Math.abs(x[i])
					: x[i];
			}
		}
		Arrays.sort(copy);
		double l = (copy.length * p) - 1;
		int i = Math.max(0,(int) Math.floor(l));
		double m = copy[i];
		return m;
	}
	
	public static float MedianSampled(float[] x,boolean fABS)
	{
		double p = 0.5;
		int len = x.length;
		int sample = 100 + (int) Math.ceil(Math.sqrt(len));
		
		float[] copy = null;
		if (sample < len / 4)
		{
			copy = new float[sample];
			for (int i = 0; i < sample; i++)
			{
				int r = (int) Math.floor(Math.random() * len);
				copy[i] = fABS
					? Math.abs(x[r])
					: x[r];
			}
		}
		else
		{
			copy = new float[len];
			for (int i = 0; i < len; i++)
			{
				copy[i] = fABS
					? Math.abs(x[i])
					: x[i];
			}
		}
		Arrays.sort(copy);
		double l = (copy.length * p) - 1;
		int i = Math.max(0,(int) Math.floor(l));
		float m = copy[i];
		return m;
	}
	
	public static float[] MedianSmooth(float[] x)
	{
		return Spectrum.MedianSmooth(x,x.length,null);
	}
	
	/** @deprecated */
	@Deprecated
	public static float[] MedianSmooth(float[] x,int len)
	{
		return Spectrum.MedianSmooth(x,len,null);
	}
	
	public static float[] MedianSmooth(float[] x,int len,float[] in)
	{
		float[] y = Spectrum.realloc(in,len);
		if (len < 3)
		{
			return x.clone();
		}
		for (int i = 1; i < len - 1; i++)
		{
			y[i] = Spectrum.Median(x[i - 1],x[i],x[i + 1]);
		}
		y[0] = y[1];
		y[len - 1] = y[len - 2];
		return y;
	}
	
	public static float[] MedianWindow(float[] x,int len,int windowSize,
		boolean fABS)
	{
		int windows = ((len - 1) / windowSize) + 1;
		float[] buckets = new float[windows];
		for (int b = 0; b < windows; b++)
		{
			int windowStart = b * windowSize;
			int windowEnd = Math.min(x.length,(b + 1) * windowSize);
			int windowLen = windowEnd - windowStart;
			buckets[b] = Spectrum.Median(x,windowStart,windowLen,fABS,null);
		}
		for (int b = 0; b < buckets.length - 1; b++)
		{
			buckets[b] = Math.min(buckets[b],buckets[b + 1]);
		}
		for (int b = buckets.length - 1; b > 0; b--)
		{
			buckets[b] = Math.min(buckets[b - 1],buckets[b]);
		}
		Spectrum.SmoothALittle(buckets);
		float[] result = Spectrum.Interpolate(buckets,len,windowSize,null);
		return result;
	}
	
	public static float[] MinimaWindow(float[] x,int len,int windowSize,
		float[] result)
	{
		result = Spectrum.realloc(result,len);
		
		float[] buckets = new float[((x.length - 1) / windowSize) + 1];
		System.arraycopy(x,0,result,0,len);
		for (int b = 0; b < buckets.length; b++)
		{
			int start = b * windowSize;
			int end = Math.min(len,(b + 1) * windowSize);
			float min = Float.MAX_VALUE;
			for (int i = start; (i < end) && (min > 0); i++)
			{
				if (result[i] < min)
				{
					min = result[i];
				}
			}
			buckets[b] = min;
		}
		for (int b = 0; b < buckets.length - 1; b++)
		{
			buckets[b] = Math.min(buckets[b],buckets[b + 1]);
		}
		for (int b = buckets.length - 1; b > 0; b--)
		{
			buckets[b] = Math.min(buckets[b - 1],buckets[b]);
		}
		Spectrum.SmoothALittle(buckets);
		result = Spectrum.Interpolate(buckets,len,windowSize,result);
		return result;
	}
	
	/* CONSIDER: Better mathematically justified ways to compute noise? */
	public static float Noise(float[] signal,int start,int end)
	{
		int len = end - start;
		float[] copy = new float[len];
		System.arraycopy(signal,start,copy,0,len);
		Arrays.sort(copy);
		int quartile = len * 3 / 4;
		return copy[quartile];
	}
	
	public static void NormalizeSquares(float[] x)
	{
		double sum = 0.0;
		for (float element : x)
		{
			sum += element * element;
		}
		double sqr = Math.sqrt(sum);
		for (int i = 0; i < x.length; i++)
		{
			x[i] /= sqr;
		}
	}
	
	public static void NormalizeSum(float[] x)
	{
		double sum = 0.0;
		for (float element : x)
		{
			sum += element;
		}
		for (int i = 0; i < x.length; i++)
		{
			x[i] /= sum;
		}
	}
	
	/**
	 * Computes points in Gaussian distribution function, over range [0,1]
	 * with mean=0
	 *
	 * @param len       number of points to compute in the distrubition (min 2)
	 * @param sigma     sqrt(variance)
	 * @param g         optional array to copy result into
	 * @return
	 */
	private static float[] NormalP(int len,double sigma,boolean even,
		boolean normalize,float[] g)
	{
		g = Spectrum.realloc(g,len);
		double variance = sigma * sigma;
		
		double delta;// = 1.0/(len-1);
		double x; // = 0;
		
		if (even)
		{
			delta = 1.0 / (len - 0.5);
			x = delta / 2.0;
		}
		else
		{
			delta = 1.0 / (len - 1.0);
			x = 0.0;
		}
		
		double s2 = -1.0 / (2.0 * variance);
		double scale = normalize
			? 1 / (sigma * Math.sqrt(2 * Math.PI))
			: 1;
		for (int i = 0; i < len; i++,x += delta)
		{
			double d = scale * Math.exp(x * x * s2);
			if (d < Float.MIN_VALUE)
			{
				break;
			}
			g[i] = (float) d;
		}
		return g;
	}
	
	public static float[] Poisson(float m)
	{
		int i = ((int) m - 5) / 10;
		i = Math.max(0,Math.min(Spectrum._poisson.length - 1,i));
		// System.err.println("**Poisson: m=" + m + ", i=" + i + ", p1=" +
		// _poisson[i][0] + ", p2=" + _poisson[i][1]);
		return Spectrum._poisson[i];
	}
	
	static int pow2(int k)
	{
		return 1 << k;
	}
	
	public static double[] realloc(double[] array,int length)
	{
		if ((null == array) || (array.length < length))
		{
			array = new double[length];
		}
		boolean debug = false;
		assert true == (debug = true);
		if (debug)
		{
			Arrays.fill(array,Double.NaN);
		}
		return array;
	}
	
	public static float[] realloc(float[] array,int length)
	{
		if ((null == array) || (array.length < length))
		{
			array = new float[length];
		}
		boolean debug = false;
		assert true == (debug = true);
		if (debug)
		{
			Arrays.fill(array,Float.NaN);
		}
		return array;
	}
	
	public static float[][] RemoveBackground(float[][] spectra)
	{
		int imzMax = spectra[0].length;
		float[][] background = new float[spectra.length][];
		for (int s = 0; s < spectra.length; s++)
		{
			float[] S = spectra[s];
			background[s] =
				Spectrum.MinimaWindow(spectra[s],spectra[s].length,72,null);
			// background[s] = Spectrum.MedianWindow(S, imzMax, 72, false);
			for (int i = 0; i < imzMax; i++)
			{
				S[i] = Math.max(0,S[i] - background[s][i]);
			}
		}
		if (spectra.length == 1)
		{
			return background;
		}
		
		float[] bg = null,row = null;
		for (int i = 0; i < imzMax; i++)
		{
			row = Spectrum.getRow(spectra,i,row);
			bg = Spectrum.MinimaWindow(row,row.length,15,bg);
			// bg = Spectrum.MedianWindow(row, row.length, 10, false);
			for (int s = 0; s < spectra.length; s++)
			{
				float b = bg[s];
				background[s][i] += b;
				spectra[s][i] = Math.max(0,spectra[s][i] - b);
			}
		}
		return background;
	}
	
	public static void Reverse(float[] x,int start,int len)
	{
		float t;
		for (int i = start,j = start + len - 1; i < j; i++,j--)
		{
			t = x[i];
			x[i] = x[j];
			x[j] = t;
		}
	}
	
	public static void Rotate(float[] x,int d)
	{
		if (d < 0)
		{
			d += x.length;
		}
		// kinda slow, but don't need another array
		Spectrum.Reverse(x,0,x.length);
		Spectrum.Reverse(x,0,d);
		Spectrum.Reverse(x,d,x.length - d);
	}
	
	public static void setRow(float[][] m,int r,float[] in)
	{
		for (int s = 0; s < m.length; s++)
		{
			m[s][r] = in[s];
		}
	}
	
	public static void SmoothALittle(double[] x)
	{
		double last = x[0];
		double cur = x[0];
		int stop = x.length - 1;
		double next;
		for (int i = 0; i < stop; i++)
		{
			next = x[i + 1];
			x[i] = (last + cur + cur + next) / 4.0F;
			last = cur;
			cur = next;
		}
		next = x[stop];
		x[stop] = (last + cur + cur + next) / 4.0F;
	}
	
	public static void SmoothALittle(float[] x)
	{
		float last = x[0];
		float cur = x[0];
		int stop = x.length - 1;
		float next;
		for (int i = 0; i < stop; i++)
		{
			next = x[i + 1];
			x[i] = (last + cur + cur + next) / 4.0F;
			last = cur;
			cur = next;
		}
		next = x[stop];
		x[stop] = (last + cur + cur + next) / 4.0F;
	}
	
	// not in Spectrum.java since this is experimental
	public static void threshold(double[][] Xin,double[] threshold)
	{
		int K = Xin.length - 1;
		
		for (int k = 0; k <= K; k++)
		{
			double[] s = Xin[k];
			int N = s.length;
			double t = threshold[k];
			if (0.0 == t)
			{
				continue;
			}
			if (Double.MAX_VALUE == t)
			{
				Arrays.fill(s,0,s.length,0.0);
				continue;
			}
			for (int i = 0; i < N; i++)
			{
				double d = s[i];
				double dm = Math.abs(d);
				s[i] = dm <= t
					? 0.0
					: d * Math.max(0,1 - Math.exp(1 - dm / t));
			}
		}
	}
	
	public static float[] UnpadToFloat(double[] y,int pad,float[] x)
	{
		int N = y.length - 2 * pad;
		if ((null == x) || (x.length != N))
		{
			x = new float[N];
		}
		for (int i = 0; i < N; i++)
		{
			x[i] = (float) y[i + pad];
		}
		return x;
	}
}
